Analyser av breareal

Importerer data

Det er fire datasett - breatlaset, 50, regioninndelingen, og omriss av Norge. Vi trenger ikke fjellmasker da vi kan anta at alle isebreer ligger i fjellet.

Breatlas

Her er det nye breatlasset med breareal fra 2018 og 2019 hentet fra Sentinell (ref Liss Marie Andreassen, NVE)

breatlas <- sf::st_read('/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/breatlas_2018_2019/Breatlas_20182019_temp20210922_lma.shp')
## Reading layer `Breatlas_20182019_temp20210922_lma' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/breatlas_2018_2019/Breatlas_20182019_temp20210922_lma.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6915 features and 27 fields (with 4 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -9874.591 ymin: 6624398 xmax: 810132 ymax: 7836994
## Projected CRS: ETRS89 / UTM zone 33N

CRS: UTM 33N ETRS89 Det er 6915 polygoner

plot(breatlas[1], main = "Breatlas 2018-2019", col = "black")

Arealet ser større ut en det er fordi det er lagt på kantlinjer.

breatlas$area <- st_area(breatlas)
par(mfrow=c(1,2))
hist(breatlas$area)
plot(breatlas$area)

De fleste polygonene er små. Og så er det noen få veldig store. Den største er 50 km2.

Kart over Norge

nor <- readRDS('../data/norway_outline.RDS')%>%
  st_as_sf()%>%
  st_transform(crs=crs(breatlas))
plot(nor$geometry, axes=T, main = "Breatlas 2018-2019")
  plot(breatlas$geometry, add=T, border = "blue")

UTM33 WGS84

N50

Her er utstrekningen til breer i perioden 1952 til 1985, basert på gamle kart.

n50 <- sf::st_read('/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/n50/cryoclim_GAO_NO_1952_1985_UTM_33N.shp')%>%
  st_transform(crs = crs(breatlas))
## Reading layer `cryoclim_GAO_NO_1952_1985_UTM_33N' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/n50/cryoclim_GAO_NO_1952_1985_UTM_33N.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 7141 features and 30 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -7278.562 ymin: 6623024 xmax: 812207.8 ymax: 7841654
## Projected CRS: WGS 84 / UTM zone 33N
plot(nor$geometry, axes=T, main = "n50 - 1952-1985")
  plot(n50$geometry, add=T, border = "red")

La oss plotte alt sammen oppå hverandre

myExt <- raster::extent(c(0, 100000, 6840000, 6900000))
myExt2 <- raster::extent(c(5000, 10000, 6880000, 6885000))

par(mfrow=c(1,3))
plot(nor$geometry, axes=T)
    plot(n50$geometry,  border = "orange", col = scales::alpha("orange", 1), add=T)
    plot(myExt, add=T, lwd=2)

plot(nor$geometry, xlim=c(0, 100000),
          ylim=c(6840000, 6900000),
          axes=T)
    plot(n50$geometry, add=T, col = "grey")
    plot(myExt2, add=T, lwd=3, col="orange")
    
plot(nor$geometry, xlim=c(5000, 10000),
          ylim=c(6880000, 6885000),
     axes=T)
    plot(n50$geometry, add=T, col = "red")
    plot(breatlas$geometry, add=T, col = "grey")

De røde områdene i kartet til høyre er hvor isen har trekt seg tilbake. Nå må vi dele kartlage inn etter regioner for så å sammenligne arealet i n50 med breatlaset. Hvordan vi gjør det med kalkulering av usikkerhet vet jeg fortsatt ikke.

Regioner

Henter shp med regionene

reg <- st_read("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Geografisk_oppdeling/regioner_2010/regNorway_wgs84 - MERGED.shp")%>%
  st_transform(crs = crs(breatlas))
## Reading layer `regNorway_wgs84 - MERGED' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Geografisk_oppdeling/regioner_2010/regNorway_wgs84 - MERGED.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4.087524 ymin: 57.75901 xmax: 31.76159 ymax: 71.38488
## Geodetic CRS:  WGS 84
plot(nor$geometry, axes=T)
  plot(reg$geometry, add=T, border = "black", 
       col = scales::alpha(c("blue", 
                             "red", 
                             "green",
                             "yellow",
                             "brown"), .2))

Før vi henter brearealet inn i tabellen for regionene må jeg prøve å fikse et problem med noen brepolygoner som krysser seg selv

table(st_is_valid(breatlas))
## 
## FALSE  TRUE 
##    68  6847
st_is_valid(breatlas, reason=T)[1:10]
##  [1] "Valid Geometry"                                  
##  [2] "Valid Geometry"                                  
##  [3] "Valid Geometry"                                  
##  [4] "Valid Geometry"                                  
##  [5] "Valid Geometry"                                  
##  [6] "Valid Geometry"                                  
##  [7] "Valid Geometry"                                  
##  [8] "Valid Geometry"                                  
##  [9] "Ring Self-intersection[100236.4931 6912507.5864]"
## [10] "Valid Geometry"
# krever lwgeom
breatlas2 <- st_make_valid(breatlas)
table(st_is_valid(breatlas2))
## 
## TRUE 
## 6915
# samme resultat:
#breatlas2 <- st_buffer(breatlas, 0.0)
#table(st_is_valid(breatlas2))
plot(nor$geometry, xlim=c(5000, 10000),
          ylim=c(6880000, 6885000),
     axes=T)
    plot(breatlas2$geometry, add=T, col = scales::alpha("blue",0.5))
    plot(breatlas$geometry, add=T, col = scales::alpha("yellow",0.5))

Dette ser ut til å ha funket fint. Denne litt kunstige inndelingen av tilgrensende polygoner er forresten bare i det nye breatlaset, ikke i n50.

breatlas <- breatlas2
rm(breatlas2)
unique(reg$region)
## [1] "Nord-Norge"      "Midt-Norge"      "Ã\u0098stlandet" "Vestlandet"     
## [5] "Sørlandet"

Fikser øæå

reg$region[reg$region=="Ã\u0098stlandet"] <- "Østlandet"
reg$region[reg$region=="Sørlandet"] <- "Sørlandet"

Breareal per region

Dagens breer

Her er en metode for å finne breareal per region.

#brealtas_reg <- st_intersection(breatlas, reg)
#saveRDS(brealtas_reg, "../data/brealtas_reg_helperfile.rds")
brealtas_reg <- readRDS("../data/brealtas_reg_helperfile.rds")

Skjekker hva som skjer med breer som ligger midt mellom to regioner

plot(brealtas_reg$geometry[brealtas_reg$region=="Midt-Norge"], 
    col = scales::alpha("red",0.5), border=NA, axes=T, 
     ylim=c(6900000, 6905000),
     xlim=c(85000, 100000))
plot(brealtas_reg$geometry[brealtas_reg$region=="Vestlandet"], 
     col = scales::alpha("green",0.5), border=NA, add=T)
  plot(nor$geometry, add=T)

Det deler polygonene ved grensen. Det er bra

brealtas_reg$area_crop <- st_area(brealtas_reg)

bretabell <- tapply(
       brealtas_reg$area_crop,
       brealtas_reg$region,
       FUN = sum)
bretabell <- bretabell/1000000 
barplot(
  bretabell,
  ylab="Breareal (km2)"
)

### Samme for N50

n50x <- st_make_valid(n50)
#n50_reg <- st_intersection(n50x, reg)
#saveRDS(n50_reg, "../data/n50_helperfile.rds")
n50_reg <- readRDS("../data/n50_helperfile.rds")
plot(n50_reg$geometry[n50_reg$region=="Midt-Norge"], 
    col = scales::alpha("red",0.5), border=NA, axes=T, 
     ylim=c(6900000, 6905000),
     xlim=c(85000, 100000))
plot(n50_reg$geometry[n50_reg$region=="Vestlandet"], 
     col = scales::alpha("green",0.5), border=NA, add=T)
  plot(nor$geometry, add=T)

n50_reg$area_crop <- st_area(n50_reg)

n50tabell <- tapply(
       n50_reg$area_crop,
       n50_reg$region,
       FUN = sum)
n50tabell <- n50tabell/1000000 


barplot(
  n50tabell,
  ylab="Breareal (km2)",
  col="grey"
)
barplot(
  bretabell,
  ylab="Breareal (km2)",
  add=T,
  col="grey30"
)

Skalerte verdier

myTbl <- tibble(region          = names(bretabell),
                indikator        = bretabell,
                referanseverdi   = n50tabell)
myTbl$skalert_indikator <- myTbl$indikator/myTbl$referanseverdi

Skalerer og klipper regionene etter norgeskartet

reg$skalert_indikator <- myTbl$skalert_indikator[match(reg$region, myTbl$region)]
reg_clipped <- st_intersection(reg, nor)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
reg_clipped$skalert_indikator_r <- round(reg_clipped$skalert_indikator, 2)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(reg_clipped) + 
  tm_polygons(col="skalert_indikator", border.col = "white")+
  tm_text("skalert_indikator_r", size=1.5)+
  tm_shape(nor)+
  tm_polygons(alpha = 0,border.col = "black")

Det er 4 måter å kalkulere den nasjonale indikatorverdien på:

  1. ta gjennomsnitt av de skalerte regionale indikatorverdiene og la alle regoinene telle likt (regional skaleringsmetode)
mean(reg_clipped$skalert_indikator)
## [1] 0.5445936
  1. summer indikatorverdien for hele landet og del på den summerte referanseverdien for hele landet (nasjonal metode)
sum(bretabell)/sum(n50tabell)
## [1] 0.7003536
  1. summere skalerte regionale indikatorverdier etter å ha gitt de ulike vekt grunnet totalarealet til regionen (arealveid metode)
reg_clipped$areal <- st_area(reg_clipped$geometry)/1000000

totalArea <- sum(reg_clipped$areal)

reg_clipped$skalert_indikator_v <- reg_clipped$skalert_indikator*
  reg_clipped$areal


sumSkalerteIndikatorer <- sum(reg_clipped$skalert_indikator_v)

sumSkalerteIndikatorer/totalArea
## 0.5884451 [1]
  1. summere skalerte regionale indikatorverdier etter å ha gitt de ulike vekt grunnet fjellareal (fjellarealveid metode)

For dette trenger jeg fjellareal per region:

(fjellareal <- readRDS("../data/fjellareal.rds"))
##       Region Fjellareal
## 1 Nord-Norge   59822.74
## 2 Midt-Norge   18294.64
## 3  Østlandet   19004.18
## 4 Vestlandet   20329.63
## 5  Sørlandet    7085.62
reg_clipped$fjellareal <- fjellareal$Fjellareal[match(reg_clipped$region, fjellareal$Region)]
reg_clipped$skalert_indikator_v2 <- reg_clipped$skalert_indikator*
  reg_clipped$fjellareal
sum(reg_clipped$skalert_indikator_v2)/sum(reg_clipped$fjellareal)
## [1] 0.619549

Alt 1 er nok ikke særlig bra. Alt 2 kan forsvares. Alt 3 er nok ikke særlig bra. Alt 4 kan forsvares og ligner veldig på hvordan NI kalkuleres.

Alt 2 sier ‘hvordan er den TOTALE ENDRINGEN i breareal i norge sett under ett’. Alt 4 sier ’hvordan er den GJENNOMSNITTLIGE endringen i breareal i i de ulike regionene i Norge.

Min intuisjon sier at alt 2 er mer riktig.

Usikkerhet

Jeg kan ikke se noen god måte å inkorporere usikkerhet inn i indikatorverdien. Det eneste jeg kan se for meg er en måte å lage en geografisk usikkerhet ved å først produsere indikatorverdier for mindre arealer (lage en raster med middels oppløning, f.eks. 10x10km2) og bootstrappe fra disse verdiene, men dette vil jeg ikke anbefale.